CL_grad_H <-
function(x, y, s, mu, A, c, inds){
	m <- dim(A)[1]
	k <- dim(A)[2]
	n <- length(inds)
	#w1bar <- x + mu*s
	#w2bar <- mu*x + s
	#wbar <- jordan_sqrt(jordan_square(w1bar) + jordan_square(w2bar) + 2*(mu^2)*jordan_identity(k))
	#Lwinv <- solve(jordan_sym(wbar))
	
	out <- matrix(0, m+k+n, m+k+n)
	out[1:m,1:k] <- -1*A
	out[(m+1):(m+k),1:k] <- CL_M(x,y,s,mu,A,c,inds)
	out[(m+1):(m+k),(k+1):(k+m)] <- -1*CL_N(x,y,s,mu,A,c,inds) %*% t(A)
	out[(m+1):(m+k),(k+m+1):(k+m+n)] <- CL_P(x,y,s,mu,A,c,inds)
	out[(m+k+1):(m+k+n),(m+k+1):(m+k+n)] <- diag(n)
	return(out)
}

